Part B

Chapter 16

  1. Using the 30-year-olds in the NHANES data set (in the NHANES package), fit a model that compares the mean height for men and women (allowing for different standard deviations). Then and answer the following questions about your model.

    1. Explain your choice of priors.
    2. Does this sample provide enough evidence to conclude that men are taller (on average)?
    3. Does this sample provide enough evidence to conclude that the standard deviation of height differs between men and women?
    Thirty <- NHANES::NHANES %>% filter(Age == 30)
    df_stats(Height ~ Gender, data = Thirty)
    ##   Gender   min    Q1 median      Q3   max     mean       sd  n missing
    ## 1 female 152.4 159.0 163.90 168.600 181.3 164.2658 6.731440 73       3
    ## 2   male 153.0 171.7 176.15 181.075 186.8 175.6744 7.082239 90       0
    gf_dens( ~ Height, color = ~ Gender, data = Thirty, binwidth = 1)
    ## Warning: Removed 3 rows containing non-finite values (stat_density).

  2. Repeat exercise 1, but compare pulse instead of height.

The previous two questions are similar. Fitting the models in JAGS is not that challenging. The interesting parts of these problems were selecting and justifying the prior and interpreting the posterior.

There is some lattitude regarding how to select a prior. Let’s suppose the likelihood is specified by

\[Y_{i|g} \sim {\sf Norm}(\mu_g, \sigma_g)\]

Then we need priors for \(\mu_g\) and \(\sigma_g\). (Technically, we need a joint prior for \(\mu_1\), \(\mu_2\), \(\sigma_1\) and \(\sigma_2\). But if the marginal priors are independent, then specifying the marginal priors determines the joint prior.)

Some important things about selecting the prior:

  • A normal prior is a fine choice for \(\mu_g\).

  • The mean of this normal prior for \(\mu_g\) should be near where we believe the mean height of people (in one group) is (prior to using our data).
    0 is not a good choice in this situation unless you use a very large \(\sigma\) because we don’t expect the mean (for eiether group) to be close to 0.

  • It is fine to use the same mean for the priors for both genders. The data will push those apart in the posterior as long as the prior isn’t too narrow.

  • The standard deviation/precision of the prior for \(\mu_g\) is not a measure of variability in the population. It is a measure of your certainty about the mean you are providing to the prior. Your explanation of how you chose the prior for the standard deviation should reflect that you understand this distinction.

  • You also need a prior for the standard deviation \(\sigma_g\). \(\sigma_g\) represents the amount of variability of \(y\) (height or pulse) in each gender group. This is the place to encorporate what you know about how large that standard deviation might be.

Note the use of separate standard deviations is not an “assumption that they are different”, but it allows them to be different and so gives us a way to test wether they are the same or different (as you were asked to do in part c).

Don’t forget to do some diagnostics. Some of you made small errors/typos that led to models that were clearly not correct but didn’t notice.

A note on the posterior: There is only one posterior (and only one prior, and one likelihood). Typically the posterior is a joint distribution, so we can talk about the marginal distributions of its components. We will often refer to these as the posterior for \(\mu\) or the posterior for \(\sigma\). It would be more explicit to refer to these as the marginal posteriors for these quantities.

When checking whether two parameters are different, look at the posterior distribution for the difference, not the amount of overlap in the two posterior distributions.

If the posterior distribution of a difference straddles 0, this does not allow you to conclude that there is no difference. If you use a ROPE and the entire HDI is contained in the ROPE, then you can conclude there is no practical difference (at least to the degree specified by the probability of your HDI).

  1. The typical lifespan of a laboratory rat that eats ad lib is approximately 700 days. When rats are placed on a restricted diet, their longevity can increase, but there is a lot of variability in lifespans across different individual rats. Restricting the diet might not only affect the typical lifespan, but restricting the diet might also affect the variance of the lifespan across rats. We consider data from R. L. Berger, Boos, and Guess (1988), as reported in Hand, Daly, Lunn, McConway, and Ostrowski (1994, data set #242), and which are available as CalvinBayes::RatLives.

    1. Run the two-group analysis on the rat longevity data using a t distribution for the response variable. Do the groups appear to differ in their central tendencies and variances? Does the value of the normality parameter suggest that the distribution of lifetimes has heavier tails than a normal distribution?

    2. Did you plot the data before doing part a? It’s a good idea to do that. Create a plot that compares the distributions of rat life for the two groups of rats in the data.

    3. Your plot should have revealed that within each group the data appear to be skewed to the left. That is, within each group, there are many rats that died relatively young, but there are fewer rats who lived especially long. We could try to implement a skewed noise distribution, or we could try to transform the data so they are approximately symmetric within each group after transformation. We will try the latter approach here. To get rid of leftward skew, we need a transformation that expands larger values more than the smaller values. We will try squaring the data. (You an use mutate() to add a new variable containing the square of the lifetimes of the rats to the original data or you can take care of this inside the list that you pass to JAGS). Do the groups appear to differ in their central tendencies and variances with this model? What does the value of the normality parameter suggest about the distribution of the transformed lifetimes?

    4. To compare the results of the two models, it is useful to back-transform to the natural scale. Give a 90% posterior HDI for the difference in mean lifetime based on each model. These should both be in units of days.

    5. Did you use ROPEs in any of your answers above? If not, go back and do so. (You will need to decide how wide the ROPE should be and if/how it changes when you apply the transformation.)

The two groups look different, but neither one looks particularly symmetric. (Other types of plots could be used here as well.)

gf_violin(DaysLive ~ Group, data = RatLives)

Nevertheless, we will proceed with this model for the moment. Let’s build a generic version of the model in JAGS that we can recycle for other purposes as we go along.

# y ~ g where y is metric and g is 1 or 2

metric_dich_model <- function() {
  for (i in 1:Nobs) {
    y[i] ~ dt(mu[g[i]], 1 / sigma[g[i]]^2, nu[g[i]])
  }
  for (g in 1:2) {
    mu[g]         ~ dnorm(mu_mean, 1 / mu_sd^2)
    sigma[g]      ~ dunif(sigma_lo, sigma_hi)
    nuMinusOne[g] ~ dexp(1 / 29.0)
    nu[g]        <- nuMinusOne[g] + 1
    log10nu[g]   <- log(nu[g]) / log(10)
  }
  delta_mu     <- mu[2] - mu[1]
  delta_sqrtmu <- sqrt(abs(mu[2])) - sqrt(abs(mu[1]))  # for sake of second attempt
  sigma_ratio  <- sigma[2] / sigma[1]
}

Now let’s fit the model.

rats_jags <-   
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "log10nu"),
    data = list(
      Nobs = nrow(RatLives),
      y = RatLives$DaysLive,
      g = as.numeric(factor(RatLives$Group)),
      mu_mean  = 700,      # approximate average life of a rat
      mu_sd    = 300,      # pretty sure mean is between 100 and 1300 (for each group)
      sigma_lo = 60 / 1000, 
      sigma_hi = 60 * 1000
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 4,
    n.thin = 2,
    DIC = FALSE
  )

Let’s take a look a what we have.

rats_jags
## Inference for Bugs model at "/var/folders/py/txwd26jx5rq83f4nn0f5fmmm0000gn/T//RtmpqhObxc/model1610a59ae5540.txt", fit using jags,
##  4 chains, each with 5000 iterations (first 1000 discarded), n.thin = 2
##  n.sims = 8000 iterations saved
##              mu.vect sd.vect    2.5%     25%      50%      75%    97.5%
## delta_mu     298.101  33.927 229.138 275.299  299.174  321.839  360.356
## log10nu[1]     0.484   0.163   0.199   0.373    0.469    0.581    0.842
## log10nu[2]     0.874   0.404   0.335   0.573    0.771    1.105    1.834
## mu[1]        706.029  10.643 684.321 698.888  706.240  713.307  726.562
## mu[2]       1004.129  32.369 937.931 982.560 1005.769 1026.805 1063.531
## sigma[1]      79.650  11.628  58.556  71.476   79.118   87.317  104.168
## sigma[2]     232.669  36.021 167.582 205.475  231.624  259.367  302.287
## sigma_ratio    2.986   0.653   1.879   2.525    2.925    3.387    4.423
##              Rhat n.eff
## delta_mu    1.001  8000
## log10nu[1]  1.001  5300
## log10nu[2]  1.002  2500
## mu[1]       1.001  8000
## mu[2]       1.001  8000
## sigma[1]    1.001  5600
## sigma[2]    1.001  8000
## sigma_ratio 1.001  8000
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
rats_mcmc <- as.mcmc(rats_jags)
mcmc_combo(rats_mcmc, regex_pars = c("delta", "ratio", "log"))
diag_mcmc(rats_mcmc)

Convergence looks good, but \(\sigma\) and \(\nu\) are being inefficiently sampled. (If we were really going to use the model, we would like run longer to boost the effective sample size.)

Suppse we are only interested in differences in average lifetime of at least 60 days, and differences in \(\sigma\) of at least 10%.

plot_post(posterior(rats_jags)$delta_mu, ROPE = c(-60, 60), xlab = "delta_mu")
## $posterior
##           ESS     mean   median     mode
## var1 2070.804 298.1006 299.1735 299.4048
## 
## $hdi
##   prob     lo       hi
## 1 0.95 232.52 362.6606
## 
## $ROPE
##    lo hi P(< ROPE) P(in ROPE) P(> ROPE)
## 1 -60 60         0          0         1
hdi(posterior(rats_jags), regex_pars = "delta_mu")
##        par     lo       hi prob
## 1 delta_mu 232.52 362.6606 0.95

plot_post(posterior(rats_jags)$sigma_ratio, ROPE = c(0.9, 1.1), xlab = "sigma_ratio")
## $posterior
##           ESS     mean   median     mode
## var1 1831.936 2.986193 2.924575 2.731553
## 
## $hdi
##   prob     lo       hi
## 1 0.95 1.7616 4.263581
## 
## $ROPE
##    lo  hi P(< ROPE) P(in ROPE) P(> ROPE)
## 1 0.9 1.1         0          0         1
hdi(posterior(rats_jags), regex_pars = "sigma_ratio")
##           par     lo       hi prob
## 1 sigma_ratio 1.7616 4.263581 0.95

hdi(posterior(rats_jags), regex_pars = "log")
##         par        lo        hi prob
## 1 log10nu.1 0.1813248 0.8099559 0.95
## 2 log10nu.2 0.2839673 1.7514411 0.95

According to this model,

  • There is a clear difference in mean longevity, with all of the posterior distribution for the difference in means is well outside a ROPE of \(\pm 2\) months (60 days).
    Indeed, the rats on the restricted diet seem to live 8 months to a year longer on average.
  • the standard deviations also appear to be different.
  • group 1 (ad lib) seems to have clearly heavier than normal tails.
  • group 2 (restricted) may also have heavier tails but the evidence is not as clear (much of the posterior credibility of \(log_10(\nu)\) is less than 1.5).

But modeling the distribution of the response as a symmetric (t) distribution doesn’t seem to align with the data. So let’s try squaring the response.

gf_violin((DaysLive^2) ~ Group, data = RatLives) 
gf_dhistogram(~ (DaysLive^2) | Group ~ ., data = RatLives, bins = 40) %>% 
  gf_fitdistr()

The squared lives do appear to be more symmetrically distributed.

Let’s a fit a new model wit the lifetimes squared. Note that the mean of a squared distribuion is not the mean of the original distribution squared:

\[E(X^2) = E(X)^2 + Var(X) \neq E(X)^2\]

rats2_jags <-   
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "delta_sqrtmu", "log10nu"),
    data = list(
      Nobs = nrow(RatLives),
      y = RatLives$DaysLive^2,
      g = as.numeric(factor(RatLives$Group)),
      mu_mean  = 700^2 + 300^2,
      mu_sd    = (700^2 + 300^2) / 3,  # takes 3 sd's to get to negative values
      sigma_lo = 500 / 1000, 
      sigma_hi = 500 * 1000
    ),
    n.iter   = 5000,
    n.burnin = 1000,
    n.chains = 4,
    n.thin   = 2,
    DIC      = FALSE
  )
summary(rats2_jags)
## fit using jags
##  4 chains, each with 5000 iterations (first 1000 discarded), n.thin = 2
##  n.sims = 8000 iterations saved
##                  mu.vect   sd.vect       2.5%        25%         50%
## delta_mu      508736.608 48178.970 415019.324 476095.302  508838.348
## delta_sqrtmu     298.580    25.430    248.680    281.440     298.647
## log10nu[1]         1.030     0.391      0.431      0.727       0.969
## log10nu[2]         1.530     0.310      0.911      1.313       1.541
## mu[1]         492920.082 15987.963 460788.700 482564.279  492914.583
## mu[2]        1001656.690 45406.261 911951.828 971722.640 1001280.963
## sigma[1]      132907.193 17554.461  99642.494 120919.095  133018.534
## sigma[2]      453425.474 26875.799 396808.997 435422.322  455529.199
## sigma_ratio        3.474     0.526      2.621      3.098       3.412
##                      75%       97.5%  Rhat n.eff
## delta_mu      540871.192  603009.022 1.002  2700
## delta_sqrtmu     315.598     347.809 1.002  2800
## log10nu[1]         1.293       1.876 1.001  8000
## log10nu[2]         1.759       2.085 1.002  2700
## mu[1]         503630.988  523900.152 1.002  2500
## mu[2]        1032322.489 1090053.088 1.002  2400
## sigma[1]      144970.909  167127.684 1.001  8000
## sigma[2]      473956.413  496266.281 1.002  2500
## sigma_ratio        3.787       4.649 1.001  8000
rats2_mcmc <- as.mcmc(rats2_jags)
plot_post(posterior(rats2_jags)$sigma_ratio, xlab = "sigma_ratio")
## $posterior
##           ESS     mean   median     mode
## var1 2237.822 3.474205 3.411965 3.345158
## 
## $hdi
##   prob       lo       hi
## 1 0.95 2.473909 4.463433
plot_post(posterior(rats2_jags)$delta_sqrtmu, xlab = "delta_sqrtmu")
## $posterior
##           ESS     mean   median   mode
## var1 6704.354 298.5803 298.6468 299.72
## 
## $hdi
##   prob       lo       hi
## 1 0.95 249.4157 348.2196

In this model we have strong evidence that the standard deviations (on the squared scale) and the means are different. The distributions of squared lifetimes appear much more normal than the untransformed lifetimes.

Let’s compare the HDIs for the difference in mean lifetimes for the two models:

HDIs <- 
  bind_rows(
    hdi(posterior(rats_jags), regex_pars = "delta_mu") %>% 
      mutate(model = "model 1"),
    hdi(posterior(rats2_jags), regex_pars = "delta_sqrtmu") %>% 
      mutate(model = "model 2")
  ) %>% 
  select(model, matches("*"))  # get model to the first column
HDIs
##     model          par       lo       hi prob
## 1 model 1     delta_mu 232.5200 362.6606 0.95
## 2 model 2 delta_sqrtmu 249.4157 348.2196 0.95
gf_dens(~delta_mu, data = posterior(rats_jags), color = ~ "model 1") %>%
  gf_dens(~delta_sqrtmu, data = posterior(rats2_jags), color = ~ "model 2") %>%
  gf_vline(xintercept = ~lo, color = ~model, data = HDIs) %>%
  gf_vline(xintercept = ~hi, color = ~model, data = HDIs) %>%
  gf_labs(x = "difference in means", color = "model") 

The centers of the posterior distributions are quite similar, but model 2 gives us a narrower HDI for the difference in the means.

  1. In the previous problem, how do the priors for the difference in mean lifetimes compare? Sample from the prior to find out. Be sure to deal appropriately with the transformation so that you are doing an apples to apples comparison.
rats_prior1 <-
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "log10nu"),
    data = list(
      Nobs = 0,
      y = RatLives$DaysLive,
      g = as.numeric(factor(RatLives$Group)),
      mu_mean  = 700,
      mu_sd    = 300,
      sigma_lo = 60 / 1000, 
      sigma_hi = 60 * 1000
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 2,
    DIC = FALSE
  )
rats_prior2 <-
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "delta_sqrtmu", "log10nu"),
    data = list(
      Nobs = 0,
      y = RatLives$DaysLive^2,
      g = as.numeric(factor(RatLives$Group)),
      mu_mean  = 700^2 + 300^2,
      mu_sd    = (700^2 + 300^2) / 3,         
      sigma_lo = 500 / 1000, 
      sigma_hi = 500 * 1000
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 2,
    DIC = FALSE
  )
gf_dens(~ delta_mu, data = posterior(rats_prior1), color = ~"model 1") %>%
gf_dens(~ delta_sqrtmu, data = posterior(rats_prior2), color = ~"model 2")

The prior for the difference in the mean lifetimes is a bit more spread out in the first model and less for the second model. If our goals was keep them similar, we should adjust one or the other. Despite this difference, both models agree that delta_mu is approximately 300. (Both priors consider 300 to be a plausible value, so it isn’t surprising that they agree in the posterior.)

  1. Shohat-Ophir et al. (2012) were interested in alcohol preferences of sexually deprived male flies. The procedure is illustrated in Figure 16.13, and was described as follows:

    One cohort, rejected-isolated, was subjected to courtship conditioning; they experienced 1-h sessions of sexual rejection by mated females, three times a day, for 4 days. …Flies in the mated-grouped cohort experienced 6-h sessions of mating with multiple receptive virgin females (ratio 1:5) for 4 days. Flies from each cohort were then tested in a two-choice preference assay, in which they voluntarily choose to consume food with or without 15% ethanol supplementation. (Shohat-Ophir et al., 2012, p. 1351, citations and figure reference removed)

    For each fly, the amount of each type of food consumed was converted to a preference ratio: the amount of ethanol-supplemented food minus the amount of regular food divided by the total of both. 3-day summary preference scores for each individual fruit fly were computed by summing the consumption of ethanol and non-ethanol across days 6–8. The amounts of food consumed and the preference ratios are in CalvinBayes::ShohatOphirKAMH2012dataReduced.

    1. How big are differences between groups relative to the uncertainty of the estimate? What do you conclude? (Answer this by computing what Kruschke calls the effect size. But note: effect size is not well defined; there are many things that go by that name. See, for example, the Wikipedia article on effect size.)

    2. Instead of focusing on the relative amounts of ethanol and regular food consumed, we might also be interested in the absolute total amount of food consumed. Run the analysis on the total consumption data, which has column name GrandTotal in the data set. What do you conclude?

flies_jags <-
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "log10nu"),
    data = list(
      Nobs = nrow(ShohatOphirKAMH2012dataReduced),
      y = ShohatOphirKAMH2012dataReduced$PreferenceIndex,
      g = as.numeric(factor(ShohatOphirKAMH2012dataReduced$Group)),
      mu_mean  = 0,      # 0 = no preference
      mu_sd    = 0.5,    # range is (-1, 1), this will sample all of that range, but prefer near 0
      sigma_lo = 0.0001, # close to 0
      sigma_hi = 1       # 1 is max possible because values in (-1, 1) 
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 2,
    DIC = FALSE
  )
## module glm loaded
## module dic loaded

(Omitting diagnositics to save space.)

Posterior distribution and HDI for delta_mu and effect size.

flies_post <- 
  posterior(flies_jags) %>% 
  mutate(efsize = delta_mu / sqrt((sigma.1^2 + sigma.2^2) / 2))
gf_dens( ~ efsize, data = flies_post)
gf_dens( ~ delta_mu, data = flies_post)
hdi(flies_post, pars = c("efsize", "delta_mu"))
##        par         lo        hi prob
## 1 delta_mu 0.08399878 0.4945389 0.95
## 2   efsize 0.25500843 2.3855492 0.95

Conclusion:

  • There appears to be a credible difference in preference (we would reject the hypothesis that there is no difference).
  • The effect size is credibly bigger than 0.2 (a benchmark often used in pyschology to indicate a “big” effect).
flies2_jags <-
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "log10nu"),
    data = list(
      Nobs = nrow(ShohatOphirKAMH2012dataReduced),
      y = ShohatOphirKAMH2012dataReduced$GrandTotal,
      g = as.numeric(factor(ShohatOphirKAMH2012dataReduced$Group)),
      mu_mean  = 150,    # rough order of maginitude of amonts eaten
      mu_sd    = 100,    # very broad
      sigma_lo = 10 / 1000,
      sigma_hi = 10 * 1000
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 2,
    DIC = FALSE
  )

Posterior distribution and HDI for delta_mu and effect size.

flies2_post <- 
  posterior(flies2_jags) %>% 
  mutate(efsize = delta_mu / sqrt((sigma.1^2 + sigma.2^2) / 2))
gf_dens( ~ efsize, data = flies2_post)
gf_dens( ~ delta_mu, data = flies2_post)
hdi(flies2_post, pars = c("efsize", "delta_mu"))
##        par          lo        hi prob
## 1 delta_mu -19.2292020 22.820278 0.95
## 2   efsize  -0.8284769  0.991917 0.95

Conclusion:

  • There is a lot of variability, but neither group clearly eats more than the other.
  • Both large and small effect sizes are credible.
  • We cannot attribute the preference difference in the previous part to differences in total consumption. (If there were a clear indication that one group consumed more, that could explain why the preference ratio was different from 0.)

Chapter 17

  1. Use Galton’s data on the men to estimate
    1. The average of height of men whose parents are 65 and 72 inches tall.
    2. The middle 50% of heights of men whose parents are 65 and 72 inches tall.
    You may use either JAGS or Stan.
set.seed(54321)
library(dplyr)
GaltonM <-
  mosaicData::Galton %>% 
  filter(sex == "M") %>%
  group_by(family) %>%
  sample_n(1) %>%
  mutate(midparent = (father + mother) / 2)
galton_model <- function() {
  for (i in 1:length(y)) {
    y[i]   ~ dt(mu[i], 1/sigma^2, nu)
    mu[i] <- alpha0 + alpha1 * (x[i] - mean(x))
  }
  sigma ~ dunif(6/100, 6 * 100)
  nuMinusOne ~ dexp(1/29)
  nu <- nuMinusOne + 1
  alpha0 ~ dnorm(0, 1/100^2)   # 100 is order of magnitude of data
  alpha1 ~ dnorm(0, 1/4^2)     # expect roughly 1-1 slope
  beta1 <- alpha1
  beta0 <- alpha0 - mean(x) * alpha1
}
library(R2jags)
library(mosaic)
galton_jags <-
  jags(
    model = galton_model,
    data = list(y = GaltonM$height, x = GaltonM$midparent),
    parameters.to.save = c("beta0", "beta1", "sigma", "nu"),
    n.iter   = 5000,
    n.burnin = 2000,
    n.chains = 4,
    n.thin   = 1
  )
Post <- 
  posterior(galton_jags) %>% 
  mutate(
    avg_ht = beta0 + beta1 * 68.5,
    ind_ht = avg_ht + rstudent_t(12000, mean = 0, sigma = sigma, df = nu))

Post %>% hdi(0.95, pars = "avg_ht")
##      par       lo       hi prob
## 1 avg_ht 70.17743 71.15686 0.95
Post %>% hdi(0.50, pars = "ind_ht")
##      par       lo       hi prob
## 1 ind_ht 69.20809 72.18985  0.5
gf_dens(~avg_ht, data = Post, color = ~ "Average") %>% 
  gf_dens(~ind_ht, data = Post, color = ~ "Individual") %>%
  gf_labs(x = "height of male with midparent = 68.5")
 
# plot_post(Post$avg_ht, hdi_prob = 0.95, xlab = "average")
# plot_post(Post$ind_ht, hdi_prob = 0.50, xlab = "individual")

  1. When centering, why did we center x but not y?

Centering on x was needed to break the correlation between the slope and intercept (which makes JAGS much more efficient). It can sometimes make the intercept easier to interpret as well.

Centering on y does not hurt anything, but the correlation has already been removed by centering on x, so centering y isn’t necessary. If we do center y, then our response essentially becomes how different is y from the average y in our data.

Note: We don’t need to center on the mean. Sometimes centering is done with convenient benchmark value. As long as the benchmark is near the mean of the values in our data, correlation will be greatly reduced. Said another way: we should pick a benchmark that is imidst the data – or design the study to collect data centered roughly on the benchmark value.

Chapter 19

  1. It turns out that larger fruit flies tend to live longer. we can take this into account by adding thorax (the length of each fruit fly’s thorax in mm) to our model. Since we are primarily interested in the association between group and longevity, thorax is referred to as a covariate. But if we were primarily interested in the relationship between longevity and size, we could use the same model and call group the covariate. Same model, different focus.

    Fit two models: Let Model 1 be the model below and let Model 2 be one just like it but without thorax. Then answer the questions that follow. We’ll stick with the default priors for now, but it would be easy enough to choose your own.

    model1_brm <-
      brm(
        bf(longevity ~ group + thorax, sigma ~ group + thorax),
        data = FruitflyReduced, 
        family = student(),
      )
    1. Does Model 1 think that longevity is different for large and small fruit flies? How do you know? How can you quantify your answer?

    2. We are primarily interested in whether the fruit flies that live with virgin females live less long than the others (those that live alone or with pregnant females). Using each model, compute a 95% HDI for the contrast that measures this. For Model 1, does it matter what value thorax has? (If yes, use thorax = 0.8 as an example.) How do the two HDI’s compare? Why?

    3. Using each model, compute a 95% HDI for the contrast that compares just the Virgin1 flies to the three types of controls (None0, Pregnant1, and Pregnant8). How do the two models compare?

    4. Now compare the values for sigma in the two models. To do this, we will need to specify both the group and the thorax size (since for different combinations of these the models estimate different values of sigma). For each of the combinations below, compute a 95% HDI for sigma in two models, one with thorax and one without. Don’t forget to convert to the natural scale. (Recall that brm() uses a log link for sigma.)

      1. group: Pregnant1; thorax: 0.7mm, 0.8mm, 0.9mm [That’s three different HDI’s for Model 1 and the same HDI 3 times for Model 2.]
      2. group: Virgin8; thorax: 0.7mm, 0.8mm, 0.9mm

      How do the results compare for the two models?

    5. Given the results above, is it important to include thorax in the model? Explain.

    6. Bonus (optional): Create a plot that shows the HDIs for sigma for every combination of model, group, and thorax measurement of 0.6, 0.7, 0.8, or 0.9 mm.

Let’s begin by fitting the models. (Be sure to use results = "hide" in your R chunk header to avoid including all of the progress indicators in your output. You may also like to use cache = TRUE so that the model isn’t refit each time you knit the document.)

model1_brm <-
  brm(
    bf(longevity ~ group + thorax, sigma ~ group + thorax),
    data = FruitflyReduced, 
    family = student(),
  )
## Compiling the C++ model
## Start sampling
model2_brm <-
  brm(
    bf(longevity ~ group, sigma ~ group),
    data = FruitflyReduced, 
    family = student(),
  )
## Compiling the C++ model
## recompiling to avoid crashing R session
## Start sampling
post1 <- posterior(model1_brm)
post2 <- posterior(model2_brm)

a The posterior for b_thorax in model 1 is well separated from 0, and the 95% HDI is

    h <- hdi(post1, pars = "b_thorax"); h
##        par       lo       hi prob
## 1 b_thorax 117.3338 158.1048 0.95

The units here is days/mm. So roughly for every 0.1 mm longer their thorax is, flies live on average between 11.7 and 15.8 days longer. (For example, the average longevity for flies with 0.8 mm thorax is roughly 11.7 to 15.8 days longer than for flies with 0.7mm).

b We can construct a contrast to measure the difference of interest. For model 1 this looks like

\[\begin{align*} \frac{\mu_{V1} + \mu_{V8}}{2} - \frac{\mu_{P1} + \mu_{P8} + \mu_{N}}{3} &= \frac{\beta_{V1} + \beta_0 + \beta_t t + \beta_{V8} + \beta_0 + \beta_t t}{2} - \frac{\beta_{P1} + \beta_0 + \beta_t t + \beta_{P8} + \beta_0 + \beta_t t + \beta_0 + \beta_t t}{3} \\ &= \frac{\beta_{V1} + \beta_{V8}}{2} - \frac{\beta_{P1} + \beta_{P8}}{3} \end{align*}\]

This will be the same for model 2 since it doesn’t have the thorax terms to begin with.

    post1 <- post1 %>%
      mutate(VvsOther = 
               (b_groupVirgin1 + b_groupVirgin8) / 2 -
               (b_groupPregnant1 + b_groupPregnant8) / 3
      )
    post2 <- post2 %>%
      mutate(VvsOther = 
               (b_groupVirgin1 + b_groupVirgin8) / 2 -
               (b_groupPregnant1 + b_groupPregnant8) / 3
      )
    bind_rows(
      hdi(post1, pars = "VvsOther") %>% mutate(model = 1),
      hdi(post2, pars = "VvsOther") %>% mutate(model = 2)
    )
##        par        lo        hi prob model
## 1 VvsOther -18.82952 -11.51699 0.95     1
## 2 VvsOther -21.49390 -10.71867 0.95     2

c

There is a lot of overlap of the two HDIs, but the HDI from model 1 is narrower (and shifted a bit lower). So model 1 provides strong evidence that the mean lifetime for the two groups differ (on average).

    post1 <- post1 %>%
      mutate(V1vsnoV = 
               (b_groupVirgin1) / 2 -
               (b_groupPregnant1 + b_groupPregnant8) / 3
      )
    post2 <- post2 %>%
      mutate(V1vsnoV = 
               (b_groupVirgin1) / 2 -
               (b_groupPregnant1 + b_groupPregnant8) / 3
      )
    bind_rows(
      hdi(post1, pars = "V1vsnoV") %>% mutate(model = 1),
      hdi(post2, pars = "V1vsnoV") %>% mutate(model = 2)
    )
##       par        lo         hi prob model
## 1 V1vsnoV -8.835859 -2.5139140 0.95     1
## 2 V1vsnoV -8.026582  0.9415611 0.95     2

d

    post1 <- post1 %>% 
      mutate(
        SigmaP1T7 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept + b_sigma_thorax * 0.7),
        SigmaP1T8 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept + b_sigma_thorax * 0.8),
        SigmaP1T9 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept + b_sigma_thorax * 0.9),
        SigmaV8T7 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept + b_sigma_thorax * 0.7),
        SigmaV8T8 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept + b_sigma_thorax * 0.8),
        SigmaV8T9 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept + b_sigma_thorax * 0.9)
        )
    post2 <- post2 %>% 
      mutate(
        SigmaP1T7 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept),
        SigmaP1T8 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept),
        SigmaP1T9 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept),
        SigmaV8T7 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept),
        SigmaV8T8 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept),
        SigmaV8T9 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept)
        )
    Sigmas <-
      bind_rows(
        hdi(post1, pars = grep("Sigma", names(post1), value = TRUE)) %>% mutate(model = 1),
        hdi(post2, pars = grep("Sigma", names(post2), value = TRUE)) %>% mutate(model = 2)
      ) %>% 
      arrange(par)
    Sigmas
##          par        lo        hi prob model
## 1  SigmaP1T7  6.264782 14.171235 0.95     1
## 2  SigmaP1T7 11.243253 20.607381 0.95     2
## 3  SigmaP1T8  8.486591 15.818850 0.95     1
## 4  SigmaP1T8 11.243253 20.607381 0.95     2
## 5  SigmaP1T9  9.589138 19.396911 0.95     1
## 6  SigmaP1T9 11.243253 20.607381 0.95     2
## 7  SigmaV8T7  3.748734  8.198663 0.95     1
## 8  SigmaV8T7  8.794834 16.404818 0.95     2
## 9  SigmaV8T8  4.999383  9.217360 0.95     1
## 10 SigmaV8T8  8.794834 16.404818 0.95     2
## 11 SigmaV8T9  5.772719 11.641738 0.95     1
## 12 SigmaV8T9  8.794834 16.404818 0.95     2
    Sigmas %>%
      gf_errorbar(lo + hi ~ par, color = ~factor(model), 
                  position = "dodge") %>%
      gf_labs(color = "model", title = "85% HPDIs for sigma")

Since thorax size explains some of the variability in lifetime, the standard deviation estimates of model 1 are smaller than those for model 2. The width of the HDI also appears smaller.

e Inclduing thorax seems to be a good idea. Longevity is positively associated with thorax size, and including thorax in the model adjusts for any possible differences in the sizes of the flies in the different groups. The smaller estimates of \(\sigma\) indicate that the model with thorax is able to give for precises estimates.

Chapter 20

  1. The CalvinBayes::Seaweed data set (adapted from [@Qian:2007]) records how quickly seaweed regenerates when in the presence of different types of grazers. Data were collected from eight different tidal areas of the Oregon coast. We want to predict the amount of seaweed from the two predictors: grazer type and tidal zone. The tidal zones are simply labeled A–H. The grazer type was more involved, with six levels: No grazers (None), small fish only (f), small and large fish (fF), limpets only (L), limpets and small fish (Lf), limpets and small fish and large fish (LfF). We would like to know the effects of the different types of grazers, and we would also like to know about the different zones.

    1. Create a plot that puts SeaweedAmt on the y-axis, Grazer on the x-axis, and uses Zone for faceting. Use gf_jitter() to avoid overplotting. Set height and width to appropriate values so the plot is still easily interpretable.

    2. Fit a model with both predictors and their interaction assuming homogeneous variances in each group. Why would it not be a good idea to fit a model with heterogenous variances?

    3. What is the effect of small fish across all the zones? Answer this question by setting up the following three contrasts: none versus small fish only; limpets only versus limpets and small fish; the average of none and limpets only versus the average of small fish only and limpets with small fish. Discuss the results.

    4. What is the effect of limpets? There are several contrasts that can address this question, but be sure to include a contrast that compares all of the locations with limpets to all of the locations without limpets.

    5. Set up a contrast to compare Zone A with Zone D. Briefly discuss the result.

    6. Does the effect of limpets depend on whether the location is in zone A or D? Use an appropriate contrast to find out.

a Looking at the data:

b Fitting a model

sw_brm <- brm(SeaweedAmt ~ Zone * Grazer, data = Seaweed)
## Compiling the C++ model
## Start sampling

Since we only have two observations in each zone/grazer group, it would be very difficult to estimate variances in each group. Since the interation model has a mean for each zone/grazer combo, if we added a separate variance for each the model would fit the data too well: If there are only two values and we know the mean and variance, we can recover the two values.

c

The biggest hassle here is creating the expression for the contrast you are interested in. I’ve written a little function to simplify this a bit, but I’ve hidden the code.

Little fish vs no little fish (no big fish in either case)

c1 <- 
  create_contrast(M[,"None"], M[, "f"], 
                  a_denom = 8, b_denom = 8, rel = "=")
c2 <- 
  create_contrast(M[,"L"], M[, "Lf"], 
                  a_denom = 8, b_denom = 8, rel = "=")
c3 <- 
  create_contrast(M[,c("L", "None")], M[, c("Lf", "f")], 
                  a_denom = 16, b_denom = 16, rel = "=")
c1
## [1] "( GrazerNone + ZoneB + GrazerNone + ZoneB:GrazerNone + ZoneC + GrazerNone + ZoneC:GrazerNone + ZoneD + GrazerNone + ZoneD:GrazerNone + ZoneE + GrazerNone + ZoneE:GrazerNone + ZoneF + GrazerNone + ZoneF:GrazerNone + ZoneG + GrazerNone + ZoneG:GrazerNone + ZoneH + GrazerNone + ZoneH:GrazerNone ) / 8 = (  + ZoneB + ZoneC + ZoneD + ZoneE + ZoneF + ZoneG + ZoneH ) / 8"
c2
## [1] "( GrazerL + ZoneB + GrazerL + ZoneB:GrazerL + ZoneC + GrazerL + ZoneC:GrazerL + ZoneD + GrazerL + ZoneD:GrazerL + ZoneE + GrazerL + ZoneE:GrazerL + ZoneF + GrazerL + ZoneF:GrazerL + ZoneG + GrazerL + ZoneG:GrazerL + ZoneH + GrazerL + ZoneH:GrazerL ) / 8 = ( GrazerLf + ZoneB + GrazerLf + ZoneB:GrazerLf + ZoneC + GrazerLf + ZoneC:GrazerLf + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneE + GrazerLf + ZoneE:GrazerLf + ZoneF + GrazerLf + ZoneF:GrazerLf + ZoneG + GrazerLf + ZoneG:GrazerLf + ZoneH + GrazerLf + ZoneH:GrazerLf ) / 8"
c3
## [1] "( GrazerL + ZoneB + GrazerL + ZoneB:GrazerL + ZoneC + GrazerL + ZoneC:GrazerL + ZoneD + GrazerL + ZoneD:GrazerL + ZoneE + GrazerL + ZoneE:GrazerL + ZoneF + GrazerL + ZoneF:GrazerL + ZoneG + GrazerL + ZoneG:GrazerL + ZoneH + GrazerL + ZoneH:GrazerL + GrazerNone + ZoneB + GrazerNone + ZoneB:GrazerNone + ZoneC + GrazerNone + ZoneC:GrazerNone + ZoneD + GrazerNone + ZoneD:GrazerNone + ZoneE + GrazerNone + ZoneE:GrazerNone + ZoneF + GrazerNone + ZoneF:GrazerNone + ZoneG + GrazerNone + ZoneG:GrazerNone + ZoneH + GrazerNone + ZoneH:GrazerNone ) / 16 = ( GrazerLf + ZoneB + GrazerLf + ZoneB:GrazerLf + ZoneC + GrazerLf + ZoneC:GrazerLf + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneE + GrazerLf + ZoneE:GrazerLf + ZoneF + GrazerLf + ZoneF:GrazerLf + ZoneG + GrazerLf + ZoneG:GrazerLf + ZoneH + GrazerLf + ZoneH:GrazerLf +  + ZoneB + ZoneC + ZoneD + ZoneE + ZoneF + ZoneG + ZoneH ) / 16"
hypothesis(sw_brm, c1) 
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerNone+Zone... = 0     9.18      3.47     2.47    16.18         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(sw_brm, c2) 
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerL+ZoneB+G... = 0     2.73      3.38    -4.07      9.4         NA
##   Post.Prob Star
## 1        NA     
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(sw_brm, c3) 
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerL+ZoneB+G... = 0     5.95      2.46     1.16    10.73         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(sw_brm, c1) %>% plot()
hypothesis(sw_brm, c2) %>% plot()
hypothesis(sw_brm, c3) %>% plot()

Limpets vs No Limpets

c4 <- 
  create_contrast(
    M[, c("None", "f", "fF")], M[, c("L", "Lf", "LfF")], 
    a_denom = 24, b_denom = 24, rel = "=")
c4
## [1] "( GrazerNone + ZoneB + GrazerNone + ZoneB:GrazerNone + ZoneC + GrazerNone + ZoneC:GrazerNone + ZoneD + GrazerNone + ZoneD:GrazerNone + ZoneE + GrazerNone + ZoneE:GrazerNone + ZoneF + GrazerNone + ZoneF:GrazerNone + ZoneG + GrazerNone + ZoneG:GrazerNone + ZoneH + GrazerNone + ZoneH:GrazerNone +  + ZoneB + ZoneC + ZoneD + ZoneE + ZoneF + ZoneG + ZoneH + GrazerfF + ZoneB + GrazerfF + ZoneB:GrazerfF + ZoneC + GrazerfF + ZoneC:GrazerfF + ZoneD + GrazerfF + ZoneD:GrazerfF + ZoneE + GrazerfF + ZoneE:GrazerfF + ZoneF + GrazerfF + ZoneF:GrazerfF + ZoneG + GrazerfF + ZoneG:GrazerfF + ZoneH + GrazerfF + ZoneH:GrazerfF ) / 24 = ( GrazerL + ZoneB + GrazerL + ZoneB:GrazerL + ZoneC + GrazerL + ZoneC:GrazerL + ZoneD + GrazerL + ZoneD:GrazerL + ZoneE + GrazerL + ZoneE:GrazerL + ZoneF + GrazerL + ZoneF:GrazerL + ZoneG + GrazerL + ZoneG:GrazerL + ZoneH + GrazerL + ZoneH:GrazerL + GrazerLf + ZoneB + GrazerLf + ZoneB:GrazerLf + ZoneC + GrazerLf + ZoneC:GrazerLf + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneE + GrazerLf + ZoneE:GrazerLf + ZoneF + GrazerLf + ZoneF:GrazerLf + ZoneG + GrazerLf + ZoneG:GrazerLf + ZoneH + GrazerLf + ZoneH:GrazerLf + GrazerLfF + ZoneB + GrazerLfF + ZoneB:GrazerLfF + ZoneC + GrazerLfF + ZoneC:GrazerLfF + ZoneD + GrazerLfF + ZoneD:GrazerLfF + ZoneE + GrazerLfF + ZoneE:GrazerLfF + ZoneF + GrazerLfF + ZoneF:GrazerLfF + ZoneG + GrazerLfF + ZoneG:GrazerLfF + ZoneH + GrazerLfF + ZoneH:GrazerLfF ) / 24"
hypothesis(sw_brm, c4) 
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerNone+Zone... = 0    28.22      1.98    24.29    32.12         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(sw_brm, c4) %>% plot()

Zone A vs Zone D

c5 <- 
  create_contrast( M["A", ], M["D", ],
    a_denom = 6, b_denom = 6, rel = "=")
c5
## [1] "(  + GrazerfF + GrazerL + GrazerLf + GrazerLfF + GrazerNone ) / 6 = ( ZoneD + ZoneD + GrazerfF + ZoneD:GrazerfF + ZoneD + GrazerL + ZoneD:GrazerL + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneD + GrazerLfF + ZoneD:GrazerLfF + ZoneD + GrazerNone + ZoneD:GrazerNone ) / 6"
hypothesis(sw_brm, c5) 
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((+GrazerfF+Graze... = 0   -45.53      3.97   -53.22   -37.77         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(sw_brm, c5) %>% plot()

c6 <- 
  create_contrast( 
    c(M["A", c("L", "Lf", "LfF")], M["D", c("None", "f", "fF")]),
    c(M["D", c("L", "Lf", "LfF")], M["A", c("None", "f", "fF")]),
    a_denom = 6, b_denom = 6, rel = "=")
c6
## [1] "( GrazerL + GrazerLf + GrazerLfF + ZoneD + GrazerNone + ZoneD:GrazerNone + ZoneD + ZoneD + GrazerfF + ZoneD:GrazerfF ) / 6 = ( ZoneD + GrazerL + ZoneD:GrazerL + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneD + GrazerLfF + ZoneD:GrazerLfF + GrazerNone +  + GrazerfF ) / 6"
hypothesis(sw_brm, c6) 
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerL+GrazerL... = 0    20.55      3.96    12.64    28.33         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(sw_brm, c6) %>% plot()

  1. This problem investigates the synthetic data set CalvinBayes::NonhomogVar.

    1. From the plot, it is pretty clear that that variance is not the same across the groups. Fit two models.
    model1_brm <- brm(Y ~ Group, data = NonhomogVar)
    model2_brm <- brm(bf(Y ~ Group, sigma ~ Group), data = NonhomogVar)
    1. What is the difference between these two models? (That is, what does sigma ~ Group do?)

    2. Describe the priors for each model. (Use prior_summary() if you are not sure.)

    3. Create the following plots:

      1. A plot showing the posterior distributions of \(\sigma\) from model 1 and each of the \(\sigma_j\)’s from model 2. (Stack multiple calls to gf_dens() on a single plot. Use color or linetype or size or some combination of these to make them distinguishable. Note: color = ~"sigmaA", etc will give you a nice legend. This works for linetype and size as well.)

      2. A plot showing the posterior distribution of \(\sigma_A - \sigma_B\) in model 2.

      3. A plot showing the posterior distribution of \(\sigma_B - \sigma_C\) in model 2.

      4. A plot showing the posterior distribution of \(\frac{\sigma_A + \sigma_D}{2} - \frac{\sigma_B + \sigma_C}{2}\) in model 2.

    4. Use each model to answer the following:

      1. Are the means of groups A and B different?
      2. Are the means of groups B and C different?

      Explain why the models agree or disagree.

    5. The original plot suggests that model 2 should be preferred over model 1. Compare the models using WAIC and LOOIC. Are these measures able to detect that model 2 is better than model 1?

model1_brm <- brm(Y ~ Group, data = NonhomogVar)
## Compiling the C++ model
## Start sampling
model2_brm <- brm(bf(Y ~ Group, sigma ~ Group), data = NonhomogVar)
## Compiling the C++ model
## Start sampling

The second of these models uses a different value of \(\sigma\) for each group. The first uses the same \(\sigma\) for all the groups.

Priors:

  • intercept (both models): \({\sf T}(3, 100, 10)\)

    df_stats(~Y, data = NonhomogVar, mean, sd)
    ##   mean_Y     sd_Y
    ## 1  100.5 6.228965

    So this is using a (order-of-magnitude rounded) mean and standard deviation of the response to determine the prior.

  • \(\log(\sigma)\) (model 1): \({\sf T}(3, 0, 10)\)
  • \(\log(\sigma_A)\) (model 2): \({\sf T}(3, 0, 10)\)
  • other linear model coefficients (including the offsets for \(\log(\sigma)\) in model 2): improper uniform prior

As expected, the model that fits separate \(\sigma\)’s for each group estimates two smaller values and to larger values. Model 1 has a single \(\sigma\) that falls between the two.

Post1 <- posterior(model1_brm)
Post2 <- posterior(model2_brm) %>%
  mutate(
    sigmaA = exp(b_sigma_Intercept),     # don't forget exp()!
    sigmaB = exp(b_sigma_Intercept + b_sigma_GroupB),
    sigmaC = exp(b_sigma_Intercept + b_sigma_GroupC),
    sigmaD = exp(b_sigma_Intercept + b_sigma_GroupD)
  ) 
Post2a <- 
  Post2 %>% 
  select(sigmaA, sigmaB, sigmaC, sigmaD) %>%
  tidyr::gather(group, sigma, sigmaA:sigmaD)

gf_dens(~sigma, data = Post1, size = 0.7) %>%
  gf_dens(~sigma, data = Post2a, color = ~group, size = 1, alpha = 0.5)

Also as expected, Model 2 finds no credible difference between the standard deviations of group B and group C, but is quite convinced that there is a difference bewteen groups A and B and between (A or D) vs (B or C).

gf_dens( ~(sigmaA - sigmaB), data = Post2, color = ~"A vs B", size = 1) %>%
gf_dens( ~(sigmaB - sigmaC), data = Post2, color = ~"B vs C", size = 1) %>%
gf_dens(~((sigmaA + sigmaD) / 2 - (sigmaB + sigmaC) / 2), color = ~ "AD vs BC",
        data = Post2, size = 1)

The HDIs tell the same story as the plots above.

bind_rows(
  hdi(Post2, pars =  "sigma: A vs B" ~ (sigmaA - sigmaB)),
  hdi(Post2, pars =  "sigma: B vs C" ~ (sigmaB - sigmaC)),
  hdi(Post2, 
      pars =  "sigma: AD vs BC" ~ ((sigmaA + sigmaD) / 2 - (sigmaB + sigmaC) / 2))
)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector
##                 par         lo       hi prob
## 1   "sigma: A vs B"  4.9898263 9.996451 0.95
## 2   "sigma: B vs C" -0.4389844 0.434171 0.95
## 3 "sigma: AD vs BC"  5.6391593 9.116402 0.95

Looking at how the two models compare the means we see that the results are quite similar when we compare group A to group B, but when we compare group B to group C, Model 2 finds compelling evidence that \(\mu_B < \mu_C\) while Model 1 does not. The reason for the difference is the Model 1 estimates each group to have a much larger standard deviation (influenced by groups A and D). The extra modeled variability within the groups hides the difference between the two groups.

bind_rows(
  hdi(Post1, pars =  "mean: B - A" ~ b_GroupB) %>% mutate(model = 1),
  hdi(Post2, pars =  "mean: B - A" ~ b_GroupB) %>% mutate(model = 2),
  hdi(Post1, pars =  "mean: B - C" ~ (b_GroupB - b_GroupC)) %>% mutate(model = 1),
  hdi(Post2, pars =  "mean: B - C" ~ (b_GroupB - b_GroupC)) %>% mutate(model = 2)
) %>% select(model, contrast = par, lo, hi, prob)
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector
##   model      contrast        lo         hi prob
## 1     1 "mean: B - A" -1.217038  5.1443755 0.95
## 2     2 "mean: B - A" -1.365880  5.1951917 0.95
## 3     1 "mean: B - C" -6.212982  0.2230876 0.95
## 4     2 "mean: B - C" -3.623713 -2.4066558 0.95

Comparing via LOO or WAIC shows that Model 2 is to be strongly preferred over Model 1. There is very little difference in the WAIC vs LOO estimates of elpd.

library(loo)
## This is loo version 2.1.0.9000.
## **NOTE: As of version 2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. Visit mc-stan.org/loo/news for details on other changes.
loo_compare(loo(model1_brm), loo(model2_brm))
##            elpd_diff se_diff
## model2_brm   0.0       0.0  
## model1_brm -64.5       8.9
loo_compare(waic(model1_brm), waic(model2_brm))
##            elpd_diff se_diff
## model2_brm   0.0       0.0  
## model1_brm -64.6       8.9
  1. Create a model like fert4_brm but use family = student(). This will add a parameter to your model. According to WAIC and LOO, which model should your prefer? (Note: the answer can be neither.)
fert4_brm <-
  brm(Yield ~ Till * Fert + (1 | Field), data = SplitPlotAgri)
## Compiling the C++ model
## Start sampling
fert4a_brm <-
  brm(Yield ~ Till * Fert + (1 | Field), data = SplitPlotAgri,
      family = student())
## Compiling the C++ model
## Start sampling

Both WAIC and LOO find the two models to be very similar in estimated out-of-sample prediction accuracy. (The differences are about the same size as the standard errors, so no strong support for one model over the other.)

loo_compare(loo(fert4_brm, reloo = TRUE), loo(fert4a_brm, reloo = TRUE))
## 8 problematic observation(s) found.
## The model will be refit 8 times.
## 
## Fitting model 1 out of 8 (leaving out observation 15)
## 
## Fitting model 2 out of 8 (leaving out observation 20)
## 
## Fitting model 3 out of 8 (leaving out observation 31)
## 
## Fitting model 4 out of 8 (leaving out observation 73)
## 
## Fitting model 5 out of 8 (leaving out observation 77)
## 
## Fitting model 6 out of 8 (leaving out observation 78)
## 
## Fitting model 7 out of 8 (leaving out observation 83)
## 
## Fitting model 8 out of 8 (leaving out observation 97)
## Start sampling
## Start sampling
## Start sampling
## Start sampling
## Start sampling
## Start sampling
## Start sampling
## Start sampling
## No problematic observations found. Returning the original 'loo' object.
##            elpd_diff se_diff
## fert4a_brm  0.0       0.0   
## fert4_brm  -1.4       0.9
loo_compare(waic(fert4_brm), waic(fert4a_brm))
##            elpd_diff se_diff
## fert4_brm   0.0       0.0   
## fert4a_brm -1.1       1.0